#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;
static doublereal c_b15 = 10.;
static doublereal c_b88 = 0.;

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dnwla_(S_fp
   op, S_fp iovect, integer *n, integer *nband, integer *nval, integer *nfig,
   integer *nperm, doublereal *val, integer *nmvec, doublereal *vec, integer
   *nblock, integer *maxop, integer *maxj, integer *nop, doublereal *p1,
   doublereal *p0, doublereal *res, doublereal *tau, doublereal *otau,
   doublereal *t, doublereal *alp, doublereal *bet, doublereal *s, doublereal
   *p2, doublereal *bound, doublereal *atemp, doublereal *vtemp, doublereal
   *d__, integer *ind, logical *small, logical *raritz, doublereal *delta,
   doublereal *eps, integer *ierr)
{
  /* System generated locals */
  integer vec_dim1, vec_offset, p0_dim1, p0_offset, p1_dim1, p1_offset, p2_dim1,
    p2_offset, t_dim1, t_offset, alp_dim1, alp_offset, bet_dim1, bet_offset,
    s_dim1, s_offset, i__1, i__2, i__3, i__4;
  doublereal d__1, d__2, d__3;

  /* Builtin functions */
  //double sqrt(doublereal), pow_dd(doublereal *, doublereal *);

  /* Local variables */
  /*extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
  integer *);*/
  static doublereal tola, temp, tolg, tmin, tmax, tarr[1];
  static logical test;
  static doublereal utol;
  //extern doublereal dnrm2_(integer *, doublereal *, integer *);
  static integer i__, j, k, l, m;
  //extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
  // integer *);
  static integer ngood, nleft;
  static doublereal anorm;
  //extern /* Subroutine */ int dmvpc_(integer *, doublereal *, integer *, 
  // integer *, doublereal *, integer *, doublereal *, doublereal *, 
  // doublereal *), dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *);
  static integer mtemp;
  static doublereal dzero[1];
  //extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
  // integer *, doublereal *, integer *);
  static integer i1;
  static doublereal pnorm, epsrt, rnorm;
  static integer ii, ng;
  //extern /* Subroutine */ int dlaeig_(integer *, integer *, integer *, 
  // integer *, doublereal *, doublereal *, integer *, doublereal *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, doublereal *), dlager_(integer *, 
  // integer *, integer *, doublereal *, doublereal *, doublereal *), 
  // dlaran_(integer *, doublereal *);
  static doublereal betmin, alpmin, betmax, alpmax;
  static integer ntheta;
  static logical enough;
  static integer number, nstart;
  //extern /* Subroutine */ int dortqr_(integer *, integer *, integer *, 
  // doublereal *, doublereal *), dvsort_(integer *, doublereal *, 
  // doublereal *, integer *, doublereal *, integer *, integer *, 
  // doublereal *);

  /* DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE */
  /* ORTHOGONALIZATION. */

  /*   NBAND  NBLOCK + 1  THE BAND WIDTH OF T. */

  /*   NVAL   THE NUMBER OF DESIRED EIGENVALUES. */

  /*   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS */
  /*     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE */
  /*     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL */
  /*     TO THE CURRENT KRYLOV SUBSPACE. */

  /*   NOP   THE NUMBER OF CALLS TO OP. */

  /*   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS. */

  /*   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. */

  /*   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. */

  /*   T   THE BAND MATRIX. */

  /*   ALP   THE CURRENT DIAGONAL BLOCK. */

  /*   BET   THE CURRENT OFF DIAGONAL BLOCK. */

  /*   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY THE BAND */
  /*      EIGENVALUE SOLVER DLAEIG. */

  /*   S   EIGENVECTORS OF T. */

  /*   SMALL   .TRUE.  IF THE SMALL EIGENVALUES ARE DESIRED. */

  /*   RARITZ  RETURNED AS  .TRUE.  IF A FINAL RAYLEIGH-RITZ PROCEDURE */
  /*     IS TO BE DONE. */

  /*   DELTA   RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE */
  /*     OF THE MATRIX.  USED IN ESTIMATING THE ACCURACY OF THE */
  /*     COMPUTED EIGENVALUES. */


  /* INTERNAL VARIABLES USED */


  /* J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT */
  /*    KRYLOV SUBSPACE. */

  /* NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS */
  /*    LIE IN THE CURRENT KRYLOV SUBSPACE). */

  /* NLEFT   THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED, */
  /*    I.E.  NLEFT = NVAL - NPERM. */

  /* NUMBER = NPERM + NGOOD. */

  /* ANORM   AN ESTIMATE OF THE NORM OF THE MATRIX. */

  /* EPS   THE RELATIVE MACHINE PRECISION. */

  /* UTOL   THE USER TOLERANCE. */

  /* TARR  AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO */
  /*       DLAEIG */

  /* DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE */
  /*       CONSISTENCY IN CALLS TO DCOPY */

  /* Parameter adjustments */
  p2_dim1 =  *n;
  p2_offset = 1+p2_dim1 * 1;
  p2 -= p2_offset;
  p0_dim1 =  *n;
  p0_offset = 1+p0_dim1 * 1;
  p0 -= p0_offset;
  p1_dim1 =  *n;
  p1_offset = 1+p1_dim1 * 1;
  p1 -= p1_offset;
  t_dim1 =  *nband;
  t_offset = 1+t_dim1 * 1;
  t -= t_offset;
  --val;
  vec_dim1 =  *nmvec;
  vec_offset = 1+vec_dim1 * 1;
  vec -= vec_offset;
  bet_dim1 =  *nblock;
  bet_offset = 1+bet_dim1 * 1;
  bet -= bet_offset;
  alp_dim1 =  *nblock;
  alp_offset = 1+alp_dim1 * 1;
  alp -= alp_offset;
  s_dim1 =  *maxj;
  s_offset = 1+s_dim1 * 1;
  s -= s_offset;
  --res;
  --tau;
  --otau;
  --bound;
  --atemp;
  --vtemp;
  --d__;
  --ind;

  /* Function Body */
  dzero[0] = 0.;
  rnorm = 0.;
  if (*nperm != 0)
  {
    /* Computing MAX */
    d__1 =  - val[1], d__2 = val[ *nperm];
    rnorm = max(d__1, d__2);
  }
  pnorm = rnorm;
  *delta = 1e31;
  epsrt = sqrt(*eps);
  nleft =  *nval -  *nperm;
  *nop = 0;
  number =  *nperm;
  *raritz = FALSE_;
  /* Computing MAX */
  d__3 = (doublereal)( - ((real)(*nfig)));
  d__1 = (doublereal)((real)(*n)) **eps, d__2 = pow_dd(&c_b15, &d__3);
  utol = max(d__1, d__2);
  j =  *maxj;

  /* ------------------------------------------------------------------ */

  /* ANY ITERATION OF THE ALGORITHM BEGINS HERE. */

  L30: i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    temp = dnrm2_(n, &p1[i__ *p1_dim1 + 1], &c__1);
    if (temp == 0.)
    {
      dlaran_(n, &p1[i__ *p1_dim1 + 1]);
    }
    /* L50: */
  }
  if (*nperm == 0)
  {
    goto L70;
  }
  i__1 =  *nperm;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    tau[i__] = 1.;
    otau[i__] = 0.;
    /* L60: */
  }
  L70: i__1 =  *n **nblock;
  dcopy_(&i__1, dzero, &c__0, &p0[p0_offset], &c__1);
  i__1 =  *nblock **nblock;
  dcopy_(&i__1, dzero, &c__0, &bet[bet_offset], &c__1);
  i__1 = j **nband;
  dcopy_(&i__1, dzero, &c__0, &t[t_offset], &c__1);
  mtemp =  *nval + 1;
  i__1 = mtemp;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    dcopy_(&j, dzero, &c__0, &s[i__ *s_dim1 + 1], &c__1);
    /* L75: */
  }
  ngood = 0;
  tmin = 1e30;
  tmax =  - 1e30;
  test = TRUE_;
  enough = FALSE_;
  betmax = 0.;
  j = 0;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. */

  L80: j +=  *nblock;

  /* THIS IS THE SELECTIVE ORTHOGONALIZATION. */

  if (number == 0)
  {
    goto L110;
  }
  i__1 = number;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (tau[i__] < epsrt)
    {
      goto L100;
    }
    test = TRUE_;
    tau[i__] = 0.;
    if (otau[i__] != 0.)
    {
      otau[i__] = 1.;
    }
    i__2 =  *nblock;
    for (k = 1; k <= i__2; ++k)
    {
      temp =  - ddot_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p1[k *p1_dim1 + 1],
        &c__1);
      daxpy_(n, &temp, &vec[i__ *vec_dim1 + 1], &c__1, &p1[k *p1_dim1 + 1],
        &c__1);

      /* THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A */
      /* NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS */
      /* TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. */

      if ((d__1 = temp * bet[k + k * bet_dim1], abs(d__1)) > (doublereal)((real)
        (*n)) *epsrt * anorm && i__ >  *nperm)
      {
        goto L380;
      }
      /* L90: */
    }
    L100: ;
  }

  /* IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. */

  L110: if (!test)
  {
    goto L160;
  }
  dortqr_(n, n, nblock, &p1[p1_offset], &alp[alp_offset]);
  test = FALSE_;
  if (j ==  *nblock)
  {
    goto L160;
  }
  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (alp[i__ + i__ * alp_dim1] > 0.)
    {
      goto L130;
    }
    m = j - (*nblock << 1) + i__;
    l =  *nblock + 1;
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      bet[i__ + k * bet_dim1] =  - bet[i__ + k * bet_dim1];
      t[l + m * t_dim1] =  - t[l + m * t_dim1];
      --l;
      ++m;
      /* L120: */
    }
    L130: ;
  }

  /* THIS IS THE LANCZOS STEP. */

  L160: (*op)(n, nblock, &p1[p1_offset], &p2[p2_offset]);
  ++(*nop);
  (*iovect)(n, nblock, &p1[p1_offset], &j, &c__0);

  /* THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      d__1 =  - bet[i__ + k * bet_dim1];
      daxpy_(n, &d__1, &p0[k *p0_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1], &c__1)
        ;
      /* L170: */
    }
    /* L180: */
  }

  /* THIS COMPUTES ALP AND P2=P2-P1*ALP. */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      ii = i__ - k + 1;
      alp[ii + k * alp_dim1] = ddot_(n, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k
        *p2_dim1 + 1], &c__1);
      d__1 =  - alp[ii + k * alp_dim1];
      daxpy_(n, &d__1, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1)
        ;
      if (k != i__)
      {
        d__1 =  - alp[ii + k * alp_dim1];
        daxpy_(n, &d__1, &p1[k *p1_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1],
          &c__1);
      }
      /* L190: */
    }
    /* L200: */
  }

  /*  REORTHOGONALIZATION OF THE SECOND BLOCK */

  if (j !=  *nblock)
  {
    goto L220;
  }
  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      temp = ddot_(n, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1);
      d__1 =  - temp;
      daxpy_(n, &d__1, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1)
        ;
      if (k != i__)
      {
        d__1 =  - temp;
        daxpy_(n, &d__1, &p1[k *p1_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1],
          &c__1);
      }
      ii = i__ - k + 1;
      alp[ii + k * alp_dim1] += temp;
      /* L210: */
    }
    /* L215: */
  }

  /* THIS ORTHONORMALIZES THE NEXT BLOCK */

  L220: dortqr_(n, n, nblock, &p2[p2_offset], &bet[bet_offset]);

  /* THIS STORES ALP AND BET IN T. */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m = j -  *nblock + i__;
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      l = k - i__ + 1;
      t[l + m * t_dim1] = alp[l + i__ * alp_dim1];
      /* L230: */
    }
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      l =  *nblock - i__ + k + 1;
      t[l + m * t_dim1] = bet[k + i__ * bet_dim1];
      /* L240: */
    }
    /* L250: */
  }

  /* THIS NEGATES T IF SMALL IS FALSE. */

  if (*small)
  {
    goto L280;
  }
  m = j -  *nblock + 1;
  i__1 = j;
  for (i__ = m; i__ <= i__1; ++i__)
  {
    i__2 = l;
    for (k = 1; k <= i__2; ++k)
    {
      t[k + i__ * t_dim1] =  - t[k + i__ * t_dim1];
      /* L260: */
    }
    /* L270: */
  }

  /* THIS SHIFTS THE LANCZOS VECTORS */

  L280: i__1 =  *nblock **n;
  dcopy_(&i__1, &p1[p1_offset], &c__1, &p0[p0_offset], &c__1);
  i__1 =  *nblock **n;
  dcopy_(&i__1, &p2[p2_offset], &c__1, &p1[p1_offset], &c__1);
  i__1 = j -  *nblock + 1;
  dlager_(&j, nband, &i__1, &t[t_offset], &tmin, &tmax);
  /* Computing MAX */
  d__1 = max(rnorm, tmax), d__2 =  - tmin;
  anorm = max(d__1, d__2);
  if (number == 0)
  {
    goto L305;
  }

  /* THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. */

  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, 
    &tmax);
  alpmin = tarr[0];
  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  dlaeig_(nblock, nblock, nblock, nblock, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, &tmax)
    ;
  alpmax = tarr[0];

  /* THIS COMPUTES ALP = BET(TRANSPOSE)*BET. */

  L305: i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      l = i__ - k + 1;
      i__3 =  *nblock - i__ + 1;
      alp[l + k * alp_dim1] = ddot_(&i__3, &bet[i__ + i__ * bet_dim1], nblock,
        &bet[k + i__ * bet_dim1], nblock);
      /* L300: */
    }
    /* L310: */
  }
  if (number == 0)
  {
    goto L330;
  }

  /* THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. */

  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  d__1 = anorm * anorm;
  dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &c_b88, 
    &d__1);
  betmin = sqrt(tarr[0]);

  /* THIS UPDATES TAU AND OTAU. */

  i__1 = number;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* Computing MAX */
    d__1 = alpmax - val[i__], d__2 = val[i__] - alpmin;
    temp = (tau[i__] *max(d__1, d__2) + otau[i__] *betmax +  *eps * anorm) /
      betmin;
    if (i__ <=  *nperm)
    {
      temp += res[i__] / betmin;
    }
    otau[i__] = tau[i__];
    tau[i__] = temp;
    /* L320: */
  }

  /* THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. */

  L330: dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  d__1 = anorm * anorm;
  dlaeig_(nblock, nblock, nblock, nblock, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &c_b88,
    &d__1);
  betmax = sqrt(tarr[0]);
  if (j <=  *nblock << 1)
  {
    goto L80;
  }

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND */
  /* LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK */
  /* IS JUSTIFIED. */

  tolg = epsrt * anorm;
  tola = utol * rnorm;
  if (*maxj - j <  *nblock ||  *nop >=  *maxop && nleft != 0)
  {
    goto L390;
  }
  goto L400;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO */
  /* SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN */
  /* ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY */
  /* IS AT 400. */

  L380: j -=  *nblock;
  *ierr =  - 8;
  L390: if (nleft == 0)
  {
    return 0;
  }
  test = TRUE_;
  L400: 
  /* Computing MIN */
  i__1 = j / 2, i__2 = nleft + 1;
  ntheta = min(i__1, i__2);
  dlaeig_(&j, nband, &c__1, &ntheta, &t[t_offset], &val[number + 1], maxj,
    &s[s_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin,  &tmax);
  dmvpc_(nblock, &bet[bet_offset], maxj, &j, &s[s_offset], &ntheta, &atemp[1],
    &vtemp[1], &d__[1]);

  /* THIS CHECKS FOR TERMINATION OF A CHECK RUN */

  if (nleft != 0 || j <  *nblock * 6)
  {
    goto L410;
  }
  if (val[number + 1] - atemp[1] > val[ *nperm] - tola)
  {
    goto L790;
  }

  /* THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T */
  /* TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED. */

  L410: if (ntheta <= nleft)
  {
    goto L470;
  }
  if (*nperm == 0)
  {
    goto L430;
  }
  m = number + nleft + 1;
  if (val[m] >= val[ *nperm])
  {
    goto L430;
  }
  --(*nperm);
  ngood = 0;
  number =  *nperm;
  ++nleft;
  goto L400;

  /* THIS UPDATES DELTA. */

  L430: m = number + nleft + 1;
  /* Computing MIN */
  d__1 =  *delta, d__2 = val[m];
  *delta = min(d__1, d__2);
  enough = TRUE_;
  if (nleft == 0)
  {
    goto L80;
  }
  ntheta = nleft;
  vtemp[ntheta + 1] = 1.;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. */

  /* THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES. */

  if (!(test || enough))
  {
    goto L470;
  }
  *delta = min(*delta, anorm);
  /* Computing MAX */
  /* Computing MAX */
  d__3 =  - val[number + 1];
  d__1 = rnorm, d__2 = max(d__3,  *delta);
  pnorm = max(d__1, d__2);
  tola = utol * pnorm;
  nstart = 0;
  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m = number + i__;
    /* Computing MIN */
    d__1 = atemp[i__] *atemp[i__] / (*delta - val[m]), d__2 = atemp[i__];
    if (min(d__1, d__2) > tola)
    {
      goto L450;
    }
    ind[i__] =  - 1;
    goto L460;

    L450: enough = FALSE_;
    if (!test)
    {
      goto L470;
    }
    ind[i__] = 1;
    ++nstart;
    L460: ;
  }

  /*  COPY VALUES OF IND INTO VTEMP */

  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    vtemp[i__] = (doublereal)((real)ind[i__]);
    /* L465: */
  }
  goto L500;

  /* THIS CHECKS FOR NEW GOOD VECTORS. */

  L470: ng = 0;
  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (vtemp[i__] > tolg)
    {
      goto L480;
    }
    ++ng;
    vtemp[i__] =  - 1.;
    goto L490;

    L480: vtemp[i__] = 1.;
    L490: ;
  }

  if (ng <= ngood)
  {
    goto L80;
  }
  nstart = ntheta - ng;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. */
  /* IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. */

  L500: test = test && !enough;
  ngood = ntheta - nstart;
  ++nstart;
  ++ntheta;

  /* THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND */
  /* EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR */
  /* FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE */
  /* EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS */
  /* PAUSE. */

  dcopy_(&ntheta, &val[number + 1], &c__1, &val[ *nperm + 1], &c__1);
  if (nstart == 0)
  {
    goto L580;
  }
  if (nstart == ntheta)
  {
    goto L530;
  }
  dvsort_(&ntheta, &vtemp[1], &atemp[1], &c__1, &val[ *nperm + 1], maxj, &j, 
    &s[s_offset]);

  /* THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING */
  /* VECTORS. */

  L530: if (!test)
  {
    nstart = 0;
  }
  if (!test)
  {
    goto L580;
  }

  /*  FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW */

  temp = atemp[1];
  i__1 = nstart;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* Computing MIN */
    d__1 = temp, d__2 = atemp[i__];
    temp = min(d__1, d__2);
    /* L535: */
  }
  m = ngood + 1;
  l = ngood + min(nstart,  *nblock);
  i__1 = l;
  for (i__ = m; i__ <= i__1; ++i__)
  {
    d__1 = temp / atemp[i__];
    dscal_(&j, &d__1, &s[i__ *s_dim1 + 1], &c__1);
    /* L540: */
  }
  m = (nstart - 1) /  *nblock;
  if (m == 0)
  {
    goto L570;
  }
  l = ngood +  *nblock;
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nblock;
    for (k = 1; k <= i__2; ++k)
    {
      ++l;
      if (l > ntheta)
      {
        goto L570;
      }
      i1 = ngood + k;
      d__1 = temp / atemp[l];
      daxpy_(&j, &d__1, &s[l *s_dim1 + 1], &c__1, &s[i1 *s_dim1 + 1],  &c__1);
      /* L550: */
    }
    /* L560: */
  }
  L570: nstart = min(nstart,  *nblock);

  /* THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. */

  L580: if (ngood == 0 || !(test || enough))
  {
    goto L600;
  }
  i__1 = ngood;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m =  *nperm + i__;
    res[m] = atemp[i__];
    /* L590: */
  }

  /* THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE */
  /* LANCZOS VECTORS. */

  L600: number =  *nperm + ngood;
  if (test || enough)
  {
    i__1 =  *n **nblock;
    dcopy_(&i__1, dzero, &c__0, &p1[p1_offset], &c__1);
  }
  if (ngood == 0)
  {
    goto L620;
  }
  m =  *nperm + 1;
  i__1 = number;
  for (i__ = m; i__ <= i__1; ++i__)
  {
    dcopy_(n, dzero, &c__0, &vec[i__ *vec_dim1 + 1], &c__1);
    /* L610: */
  }
  L620: i__1 = j;
  i__2 =  *nblock;
  for (i__ =  *nblock; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    (*iovect)(n, nblock, &p2[p2_offset], &i__, &c__1);
    i__3 =  *nblock;
    for (k = 1; k <= i__3; ++k)
    {
      m = i__ -  *nblock + k;
      if (nstart == 0)
      {
        goto L640;
      }
      i__4 = nstart;
      for (l = 1; l <= i__4; ++l)
      {
        i1 = ngood + l;
        daxpy_(n, &s[m + i1 * s_dim1], &p2[k *p2_dim1 + 1], &c__1, &p1[l
          *p1_dim1 + 1], &c__1);
        /* L630: */
      }
      L640: if (ngood == 0)
      {
        goto L660;
      }
      i__4 = ngood;
      for (l = 1; l <= i__4; ++l)
      {
        i1 = l +  *nperm;
        daxpy_(n, &s[m + l * s_dim1], &p2[k *p2_dim1 + 1], &c__1, &vec[i1
          *vec_dim1 + 1], &c__1);
        /* L650: */
      }
      L660: ;
    }
    /* L670: */
  }
  if (test || enough)
  {
    goto L690;
  }

  /* THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE */
  /* TAU RECURRENCE. */

  m =  *nperm + 1;
  i__2 = number;
  for (i__ = m; i__ <= i__2; ++i__)
  {
    temp = 1. / dnrm2_(n, &vec[i__ *vec_dim1 + 1], &c__1);
    dscal_(n, &temp, &vec[i__ *vec_dim1 + 1], &c__1);
    tau[i__] = 1.;
    otau[i__] = 1.;
    /* L680: */
  }

  /*  SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG */

  dcopy_(&ntheta, &val[ *nperm + 1], &c__1, &vtemp[1], &c__1);
  dvsort_(&ntheta, &vtemp[1], &atemp[1], &c__0, tarr, maxj, &j, &s[s_offset]);
  goto L80;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE */
  /* PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING */
  /* THE PERMANENT VECTORS. */

  L690: if (ngood == 0 &&  *nop >=  *maxop)
  {
    goto L810;
  }
  if (ngood == 0)
  {
    goto L30;
  }

  /* THIS ORTHONORMALIZES THE VECTORS */

  i__2 =  *nperm + ngood;
  dortqr_(nmvec, n, &i__2, &vec[vec_offset], &s[s_offset]);

  /* THIS SORTS THE VALUES AND VECTORS. */

  if (*nperm != 0)
  {
    i__2 =  *nperm + ngood;
    dvsort_(&i__2, &val[1], &res[1], &c__0, &temp, nmvec, n, &vec[vec_offset]);
  }
  *nperm += ngood;
  nleft -= ngood;
  /* Computing MAX */
  d__1 =  - val[1], d__2 = val[ *nperm];
  rnorm = max(d__1, d__2);

  /* THIS DECIDES WHERE TO GO NEXT. */

  if (*nop >=  *maxop && nleft != 0)
  {
    goto L810;
  }
  if (nleft != 0)
  {
    goto L30;
  }
  if (val[ *nval] - val[1] < tola)
  {
    goto L790;
  }

  /* THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED */
  /* TO LOOK FOR UNDISCLOSED MULTIPLICITIES. */

  m =  *nperm -  *nblock + 1;
  if (m <= 0)
  {
    return 0;
  }
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    l = i__ +  *nblock - 1;
    if (val[l] - val[i__] < tola)
    {
      goto L30;
    }
    /* L780: */
  }

  /* THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ */
  /* PROCEDURE IS NEEDED. */

  L790: m =  *nperm -  *nblock;
  if (m <= 0)
  {
    return 0;
  }
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    l = i__ +  *nblock;
    if (val[l] - val[i__] >= tola)
    {
      goto L800;
    }
    *raritz = TRUE_;
    return 0;
    L800: ;
  }

  return 0;

  /* ------------------------------------------------------------------ */

  /* THIS REPORTS THAT MAXOP WAS EXCEEDED. */

  L810: *ierr =  - 2;
  goto L790;

} /* dnwla_ */
